################################################################################
# Scripts for a tonnage model
################################################################################

##
## Preparatory steps
##

# Load the package
library(MapMark4)

# Set the working directory
oldWd <- setwd("D:\\Datap-Qmra\\MapMark4\\Temp\\TestCalculationsTm")

write.csv(ExampleTm, file = "myTm.csv", row.names = FALSE)
write.csv(ExampleDepEst2, file = "myDepEst.csv", row.names = FALSE)

##
## Probability calculations
##

# Generate the meta data
oMeta <- Metadata("PT002", "myTm.csv", "myDepEst.csv",
                  seed = 777, tractName = "Many Lucky Strikes",
                  depositModel = "Computer generated tm with uranium",
                  personName = "John Doe",
                  otherInfo = "This is an example assessment.")
summary(oMeta)

# Generate the pdf that represents the contained metal tonnages and check it
tm <- read.csv(getModelFilename(oMeta))
oTonPdf <- TonnagePdf(tm, seed = getSeed(oMeta))
plot(oTonPdf)
summary(oTonPdf)

# Generate the pmf that represents the number of undiscovered deposits
# and check it
depEst <- read.csv(getDepEstFilename(oMeta))
oPmf <- NDepositsPmf("NegBinomial", list(nDepEst = depEst))
plot(oPmf, areLinesAdded = FALSE)
summary(oPmf)

# Simulate undiscovered deposits in the permissive tract
oSimulation <- Simulation(oPmf, oTonPdf, seed = getSeed(oMeta))

# Generate the pdf that represents the total contained metal tonnages for all
# simulated undiscovered deposits in the permissive tract and check it
oTotalTonPdf <- TotalTonnagePdf(oSimulation, oPmf, oTonPdf)
plot(oTotalTonPdf)
plotMarginals(oTotalTonPdf)
summary(oTotalTonPdf)

##
## Archive of the calculation results
##

# Summaries of the probability calculations

filename <- paste(getTractId(oMeta), "_summary.txt", sep = "")
oldWidth <- options(width = 120)
sink(filename)
summary(oMeta)
summary(oTonPdf)
summary(oPmf)
summary(oTotalTonPdf)
sink()
options(width = oldWidth$width)

# Simulated undiscovered deposits, which are needed for
# subsequent assessment calculations that are not part of MapMark4
filename <- paste(getTractId(oMeta), "_Simulation.csv", sep = "")
print(oSimulation, filename)

# Instantiated classes
filename <- paste(getTractId(oMeta), ".dat", sep = "")
save(oMeta, oTonPdf, oPmf, oSimulation, oTotalTonPdf, file = filename)

# Store the plots in PNG format
png(filename = paste(getTractId(oMeta), "_TonPdf.png", sep = ""), height = 240)
plot(oTonPdf)
dev.off()

png(filename = paste(getTractId(oMeta), "_Pmf.png", sep = ""), height = 240)
plot(oPmf, areLinesAdded = FALSE)
dev.off()

png(filename = paste(getTractId(oMeta), "_TotalTonPdf.png", sep = ""), height = 240)
plot(oTotalTonPdf)
dev.off()

png(filename = paste(getTractId(oMeta), "_TotalTonPdf_Marg.png", sep = ""))
plotMarginals(oTotalTonPdf)
dev.off()


# Store the plots in postscript format
postscript(paste(getTractId(oMeta), "_TonPdf.ps", sep = ""), height = 3.5)
plot(oTonPdf)
dev.off()

postscript(paste(getTractId(oMeta), "_Pmf.ps", sep = ""), height = 3.5)
plot(oPmf)
dev.off()

postscript(paste(getTractId(oMeta), "_TotalTonPdf.ps", sep = ""), height = 3.5)
plot(oTotalTonPdf)
dev.off()

postscript(paste(getTractId(oMeta), "_TotalTonPdf_Marg.ps", sep = ""))
plotMarginals(oTotalTonPdf)
dev.off()

setwd(oldWd)